BMX: Biological modelling and interface exchange

High performance computing has a great potential to provide a range of significant benefits for investigating biological systems. These systems often present large modelling problems with many coupled subsystems, such as when studying colonies of bacteria cells. The aim to understand cell colonies has generated substantial interest as they can have strong economic and societal impacts through their roles in in industrial bioreactors and complex community structures, called biofilms, found in clinical settings. Investigating these communities through realistic models can rapidly exceed the capabilities of current serial software. Here, we introduce BMX, a software system developed for the high performance modelling of large cell communities by utilising GPU acceleration. BMX builds upon the AMRex adaptive mesh refinement package to efficiently model cell colony formation under realistic laboratory conditions. Using simple test scenarios with varying nutrient availability, we show that BMX is capable of correctly reproducing observed behavior of bacterial colonies on realistic time scales demonstrating a potential application of high performance computing to colony modelling. The open source software is available from the zenodo repository https://doi.org/10.5281/zenodo.8084270 under the BSD-2-Clause licence.

. A slight discontinuity may be seen at 200000 seconds due to a checkpoint restart. b) Log-log plot of speedup curve for system shown in Figure 1. The curve for ideal speedup is included for reference. lower performance than using the optimum number of GPUs.
It is also instructive to compare the performance of the two main components of the simulation, the computation on the grid and the computation over particles for GPUs and CPUs. Figure S.2-b shows the amount of time per time step spent in the grid and particle updates. For calculations done on the CPU, the two components are roughly comparable, the particle updates taking less than twice the time spent in the grid updates. On GPUs, however, the particle updates are nearly an order of magnitude faster than the grid updates. Although the chemistry model in the particles is extremely simple, these results suggest that the BMX framework will be able to support much more complicated chemistry inside the particles, provided it can be kept on the GPUs.

S.2 Transport Model
As discussed in the main article, the transport of nutrients through the nutrient medium is governed by the diffusion equation where the symbols were defined previously in the main article. Based on earlier work by Weissberg 3 , Khirevich et al. tried mapping their results for diffusion in packed beds using the following form for the tortuosity where p is an adjustable parameter. Weissberg originally proposed a value of 0.5 for p. Depending on the details of the packing algorithm used to create the packed bed, Khirevich et al. found values of p that varied from 0.43 to 0.5. They also found fairly high values of R 2 for their fits. From this they concluded that this form of the tortuosity was insufficient for their purposes (modeling chromatography columns). However, given the uncertainties in modeling biological systems, we feel that this form is accurate enough to model extracellular diffusion. We used a value of 0.45 for p. A problem with this formula is that it is generally possible for the volume fraction of fluid in the grid cell to go to zero, resulting in negative or extremely small diffusion coefficients. This is undesirable, since it completely shuts off transport of nutrients into and out of the particles in this grid cell from the surrounding fluid. The particle stops growing, even under circumstances where we would expect that nutrients could still get to the particle through mechanisms such as diffusion through a thin fluid film on the cell surface or diffusion through other cells. To prevent unphysical low diffusion, D e f f ,i is not allowed to drop below 0.1D 0 i .

S.3 Particle Model
The particles in the simulations described here each possess a simple "ABC" metabolic model. The chemical species in the metabolism are related to each other via the following reactions: As described previously, to distinguish between concentrations of these components and absolute amounts of the component, we use brackets [ ]. In this notation, [A out ] is the concentration of A in a grid cell and A out is the total amount (mass) of A in a grid cell. The concentrations of these species are governed by the set of equations These reactions also act as sources and sinks for the diffusion equations governing transport of A out and C out in the agar growth medium. These coupled rate equations could be handled by a time integration package but we chose to implement a simple integration scheme for two reasons. The first is that given the data layout in AMReX and their model for interacting to GPUs it would be difficult or impossible to incorporate a library-base time integration package into the code. The second reason is that it is not clear that the concentration positivity constraints, discussed in the next section, for concentrations both in the cell and in the fluid could be easily enforced using a standard integrator.
The mass flux equations into and out of the particle can be examined in more detail. The particle membrane is permeable only to the chemicals A and C and we define coupled coefficients as mass-transfer coefficients for transport of A and B, respectively. The fluxes of material into the cell, f in , and out of the cell, f out , can be written as In a time increment ∆t, the chemical mass increments ∆A in and ∆C in going into the particle and ∆A out and ∆C out going out of the particle are Note that (A in , A out ,C in ,C out ) are the absolute amounts of the chemicals occupying a particle volume V particle and fluid volume in a grid cell V grid , respectively. We are ignoring possible changes in the particle and grid cell volumes during this time interval. The incremental changes in the total amount of material in the particle and the grid cells can be written as Dividing the first equation by ∆tV cell and the second by ∆tV grid and taking the limit ∆t → 0 gives Using these expressions for the change in concentration due to transport across the cell membrane, the rate equations governing this system become From here on in, we drop the primes on k ′ 1 , k ′ −1 , k ′ 3 and k ′ −3 . Numerically, the rate equations are currently being handled using the following scheme: • There is an equilibrium between the concentration of species in the fluid and the particle over a time increment ∆t/2 • The internal particle reactions take place over a time increment ∆t using the updated internal concentrations from the previous step. The cell volume is also updated during this step • Using the updated internal concentrations from the reaction step, the concentrations between the inside and outside fluids are re-equilibrated over a time interval ∆t/2 The first step is implemented by calculating the amount of material that flows between the particle and the grid cell in which it is located. The total amount of material that flows into the cell in the interval ∆t is Based on this increment, the concentrations inside and outside the particle are adjusted using the equations The second step adjusts the concentration of reactants inside the particle based on internal chemical reactions. Using a simple Euler scheme, the adjusted values of the concentrations and the particle volume after the time increment ∆t are After adjusting the volumes, the concentrations need to be modified slightly by multiplying by the ratio of V particle and V ′ particle The concentrations inside and outside the particle are then adjusted one more time using the equations in the first step. The increments in A and C between the exterior fluid and the particle interior from the first and last step are added together to give a total increment between the fluid and the particle interior. At the end of the step, the increment is divided by the fluid volume, V grid , in the grid cell to get the change in concentration in the grid cell hosting the particle.

S.4 Maintaining Non-negative Concentrations
Chemical concentrations are positive quantities. Therefore an important part of successfully running a simulation is to guarantee that at no time do the concentrations of any chemical species fall below zero, in either the particle or in the grid cell. This can occur because the finite value of the time increment can cause values to overshoot their true values by a small amount, resulting in negative concentrations after exchanges. The potential for error is amplified when considering grid volume that contain multiple particles, each simultaneously adding or subtracting from the grid concentrations. This can be rectified by using a smaller time step, but this is also undesirable since it results in longer simulation times with relatively little change in physical behavior. A better approach is to modify the algorithm to correct for negative concentrations when they occur.
The BMX strategy for maintaining positive concentrations consists of two parts. The first is to divide the available fluid in each grid cell evenly between each of the particles located in the grid cell and the second is to check if any of the components in the reaction schemes inside the particles goes negative. Apportioning the fluid volume evenly between the particles is designed to handle the issue that the chemistry of multiple particles can potentially interact by exchanging material through the fluid that they all share in common. However, a chemical integration scheme that updates particles individually will not account for this and while each individual particle may be able to maintain non-negative concentrations, multiple particles acting independently of each other in a single grid cell could result in a negative fluid concentration. If the fluid is divided up between the each particle in the grid cell, then this cannot happen if the chemistry of the individual particles has been implemented to guarantee that the internal concentration and fluid concentrations all remain positive.
Concentrations inside the particle or in the fluid can go negative because a finite time step causes the system to overshoot a limit on the available reactant. If this occurs, then the change in that component is modified so that the component only goes to zero. This means that the change in any component that is coupled to the changed component must also be adjusted. For example, for the reaction we might calculate a change in A over one time step, ∆A, such that A goes negative. This can be corrected by decreasing the change in A to a value ∆A ′ so that A goes to zero instead. However, this means that the changes in B and C must also be modified to reflect the smaller change in A.

S.4.1 Parameters
In this article we have viewed the particles as representing bacteria cells utilising a simple "ABC" metabolic model. Even with the simple "ABC" model presented here, this model has a substantial number of parameters that need to be set in order for the simulation to produce meaningful results. For the purposes of this study, the parameters were chosen to give results comparable to the behavior of a typical bacterial colony with the qualitative results being compared to literature and experimental observations. The parameters used were sufficient to demonstrate that the growth model can be handled by the existing numerical framework and integration schemes introduced within BMX. We expect that a more realistic metabolism model would be characterized by many more reactions and parameters with increased simulation complexity, but from a numerical point of view would not be qualitatively different from the model presented here.
The volume of a bacteria cell that is 1 µm in diameter is 5.3×10 −13 cm 3 so the volume threshold for splitting into 2 cells, V crit , was set at approximated twice that value, 1×10 −12 cm 3 (the radius at splitting is about 6.2 µm). Given that the dimension of individual cells are in the range of 1 µm in diameter, a minimum grid cell size was also set to 1 µm. The remaining parameters were chosen so that particle division occurred on a time scale of 30 minutes to an hour. The diffusion coefficient for a nutrient molecule such as glucose is approximately 600 cm 2 /s and this was used as the diffusion coefficient for both D A and D C in the extracellular medium. The initial concentration of A was chosen to be 20 mMolar. The concentration of C was chosen to have a non-zero value at the beginning of the simulation and was arbitrarily set to 2 mMolar. The component B is only found inside the bacteria cell particles and its concentration in the extracellular medium is zero.
The rate constants were chosen through a combination of analysis and empirical adjustment to give the desired behavior. Values were adjusted to give a timescale of approximately 30 minutes to an hour for cell division. For a realistic model, these parameters would be known, or mostly known, a priori and would not be subjected to additional adjustment. The complete set of parameters is summarized in

S.5 Installation and Usage Instructions
BMX has been developed to provide a portable high-performance agent-based simulation suite for large scale cell community modelling. The portable nature of BMX means that the software can be installed and run on a range of different hardware architectures with minimal user intervention. One of the advancing capabilities of BMX is the GPU acceleration which is shown in Section 1 to greatly improve the computation speed for large cell community simulations. A GPU is not however a requirement to run BMX and good simulation speeds can be achieved through CPU computation alone. Therefore, BMX can be readily used across a range of computing systems and as such provides an accessible suite for large scale simulations. Adaptability was a design aim for BMX allowing a range of biological problems to be considered. A user's simulation may be implemented by modifying the input files and by editing the chemical evolution routines of the c++ source code. The input files are used to determine both the solver/simulation settings, to provide the rate parameters for chemical reactions, and input the simulation initial conditions. The simulation outputs are provided as standardized .vtk files for rapid post-processing by a user's preferred visualization software.